function hfile = vmp_ImportSPMMaps(hfile, maps, opts)
% VMP::ImportSPMMaps  - import SPM maps from Analzye/NIftI files
%
% FORMAT:       vmp.ImportSPMMaps(maps [, opts])
%
% Input fields:
%
%       maps        filenames of spmT_* / spmF_* maps
%       opts        1x1 struct with optional fields
%        .interp    method, one of 'nearest', {'linear'}, 'cubic'
%        .maptype   1x1 or Mx1 map type, 'F', 'r', 't' (default from name)
%        .thresh    1x2 or Mx2 lower and upper thresholds (default: p=0.05)
%
% No output fields.

% Version:  v0.7f
% Build:    9030323
% Date:     Mar-03 2009, 11:06 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmp') || ...
    isempty(maps) || ...
   (~ischar(maps) && ...
    ~iscell(maps))
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
res = bc.Resolution;
if isempty(bc.Map)
    error( ...
        'BVQXfile:BadArgument', ...
        'The VMP must not be empty (for default settings).' ...
    );
end
if nargin < 3 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ischar(maps)
    if size(maps) == 1
        maps = {maps(:)'};
    else
        maps = cellstr(maps);
    end
end
if ~isfield(opts, 'interp') || ...
   ~ischar(opts.interp) || ...
   ~any(strcmpi(opts.interp(:)', {'cubic', 'linear', 'nearest'}))
    opts.interp = 'linear';
else
    opts.interp = lower(opts.interp(:)');
end
if ~isfield(opts, 'maptype') || ...
   ~ischar(opts.maptype) || ...
   (numel(opts.maptype) ~= 1 && ...
    numel(opts.maptype) ~= numel(maps))
    opts.maptype = 'a';
else
    opts.maptype = lower(opts.maptype(:)');
end
if numel(opts.maptype) == 1
    opts.maptype(2:numel(maps)) = opts.maptype(1);
end
if ~isfield(opts, 'thresh') || ...
   (~isa(opts.thresh, 'double') && ...
    ~isa(opts.thresh, 'single')) || ...
    size(opts.thresh, 2) ~= 2 || ...
   (size(opts.thresh, 1) ~= 1 && ...
    size(opts.thresh, 1) ~= numel(maps)) || ...
    any(isinf(opts.thresh(:)) | isnan(opts.thresh(:)))
    opts.thresh = [0.05, 0.001];
else
    opts.thresh = abs(opts.thresh(:, 1:2, 1));
    opts.thresh = opts.thresh(:, 1:2);
end
if size(opts.thresh, 1) == 1
    opts.thresh = repmat(opts.thresh, [numel(maps), 1]);
end
for mc = numel(maps):-1:1
    if ~ischar(maps{mc}) || ...
        isempty(maps{mc}) || ...
        exist(maps{mc}(:)', 'file') ~= 2
        maps(mc) = [];
        opts.maptype(mc) = [];
        opts.thresh(mc, :) = [];
        continue;
    end
    maps{mc} = maps{mc}(:)';
    if numel(maps{mc}) > 4 && ...
        strcmpi(maps{mc}(end-3:end), '.img')
        if exist([maps{mc}(1:end-3) 'hdr'], 'file') == 2
            maps{mc} = [maps{mc}(1:end-3) 'hdr'];
        elseif exist([maps{mc}(1:end-3) 'HDR'], 'file') == 2
            maps{mc} = [maps{mc}(1:end-3) 'HDR'];
        else
            opts.maptype(mc) = [];
            opts.thresh(mc, :) = [];
            maps(mc) = [];
        end
    end
end
if isempty(maps)
    error( ...
        'BVQXfile:BadArgument', ...
        'Map file(s) not found.' ...
    );
end
nobj = numel(maps);
mobj = cell(1, nobj);
for mc = 1:nobj
    try
        mobj{mc} = BVQXfile(maps{mc});
        if numel(mobj{mc}) ~= 1 || ...
           ~isBVQXfile(mobj{mc}, 'hdr')
            error('BAD_HDR_FILE');
        end
        tbc = bvqxfile_getcont(mobj{mc}.L);
        if isempty(regexpi(tbc.DataHist.Description, ...
                '^spm (beta|contrast) - (\d+)\:\s*(.*)$')) && ...
            isempty(regexpi(tbc.DataHist.Description, ...
                '^spm\{([a-z]+)_\[([0-9][0-9\.]+)(\,[0-9][0-9\.]+)?\]\}\s*\-?\s*(.*)$'))
            if opts.maptype(mc) == 'a'
                warning( ...
                    'BVQXfile:AutoDetectFailed', ...
                    'Error detecting SPM image type, setting to beta.' ...
                );
                tbc.DataHist.Description = 'spm beta - 0: imported map';
            end
        end
    catch
        clearbvqxobjects(mobj);
        error( ...
            'BVQXfile:BadArgument', ...
            'Invalid map file (not an SPM result map): %s', ...
            maps{mc} ...
        );
    end
end

% build coordinate indices for sampling
bb = struct('BBox', [bc.XStart, bc.YStart, bc.ZStart; bc.XEnd, bc.YEnd, bc.ZEnd], 'XYZRes', res);
if res > 1 || ...
    bc.NativeResolutionFile
    bb.BBox(2, :) = bb.BBox(2, :) - 1;
else
    bb.BBox(2, :) = bb.BBox(2, :) + 0.5;
end

% add maps to structure once !
bc.Map = bc.Map(:);
onmaps = numel(bc.Map);
bc.Map(onmaps+1:onmaps+nobj) = bc.Map(onmaps);

% FDR thresholds
fdrv = [0.1, 0.05, 0.04, 0.03, 0.02, 0.01, 0.005, 0.001]';

% iterate over map objects
for oc = 1:nobj
    
    % get target position
    tobj = onmaps + oc;
    
    % get content of map
    sobj = bvqxfile_getscont(mobj{oc}.L);
    cobj = sobj.C;
    if istransio(cobj.VoxelData)
        cobj.VoxelData = resolve(cobj.VoxelData);
    end
	cobj.VoxelData = single(cobj.VoxelData);
	cobj.VoxelData(isnan(cobj.VoxelData)) = single(0);
    bvqxfile_setcont(mobj{oc}.L, cobj);
    bfv = sum(cobj.VoxelData(:) ~= 0);
    
    % parse description
    [sfn{1:2}] = fileparts(sobj.F);
    if ~isempty(sfn{1})
        [pfn{1:2}] = fileparts(sfn{1});
        sfn{2} = [pfn{2} '/' sfn{2}];
    end
    desc = cobj.DataHist.Description;
    [iti{1:3}] = regexpi(desc, '^spm (beta|contrast) - (\d+)\:\s*(.*)$');
    [dfi{1:3}] = regexpi(desc, ...
        '^spm\{([a-z]+)_\[([0-9][0-9\.]+)(\,[0-9][0-9\.]+)?\]\}\s*\-?\s*(.*)$');
    
    % maptype
    if ~any('bcfrtz' == opts.maptype(oc))
        if ~isempty(iti{3})
            opts.maptype(oc) = lower(desc(iti{3}{1}(1)));
            mapname = sprintf('%s - %s %s: %s', sfn{2}, ...
                desc(iti{3}{1}(1, 1):iti{3}{1}(1, 2)), ...
                desc(iti{3}{1}(2, 1):iti{3}{1}(2, 2)), ...
                desc(iti{3}{1}(3, 1):iti{3}{1}(3, 2)));
        elseif ~isempty(dfi{3})
            opts.maptype(oc) = lower(desc(dfi{3}{1}(1)));
            mapname = sprintf('%s: %s', sfn{2}, ...
                desc(dfi{3}{1}(4, 1):dfi{3}{1}(4, 2)));
        else
            mapname = sprintf('%s: %s', sfn{2}, desc);
        end
    else
        mapname = sprintf('%s: %s', sfn{2}, desc);
    end
    if ~any('bcfrtz' == opts.maptype(oc))
        opts.maptype(oc) = 'b';
    end
    
    % fill in map values
    bc.Map(tobj).DF2 = 0;
    switch (opts.maptype(oc))
        case {'b', 'c'}
            bc.Map(tobj).Type = 11;
        case {'f'}
            bc.Map(tobj).Type = 4;
            if ~isempty(dfi) && ...
               ~isempty(dfi{3}) && ...
            	dfi{3}{1}(3, 2) >= dfi{3}{1}(3, 1)
                bc.Map(tobj).DF2 = ...
                    str2double(desc(dfi{3}{1}(3, 1):dfi{3}{1}(3, 2)));
            else
                bc.Map(tobj).DF2 = 1;
            end
        case {'r'}
            bc.Map(tobj).Type = 2;
        case {'t'}
            bc.Map(tobj).Type = 1;
        case {'z'}
            bc.Map(tobj).Type = 12;
    end
    if any('frt' == opts.maptype(oc))
        if ~isempty(dfi{3})
            bc.Map(tobj).DF1 = str2double(desc(dfi{3}{1}(2, 1):dfi{3}{1}(2, 2)));
        else
            bc.Map(tobj).DF1 = 120;
        end
        bc.Map(tobj).FDRThresholds = [fdrv, ...
            applyfdr(double(cobj.VoxelData(:)), opts.maptype(oc), fdrv, ...
            bc.Map(tobj).DF1, bc.Map(tobj).DF2, true)];
    elseif (opts.maptype(oc) == 'z');
        bc.Map(tobj).DF1 = 5e3;
        bc.Map(tobj).FDRThresholds = [fdrv, ...
            applyfdr(double(cobj.VoxelData(:)), 't', fdrv, ...
            bc.Map(tobj).DF1, bc.Map(tobj).DF2, true)];
    else
        bc.Map(tobj).DF1 = 1;
        bc.Map(tobj).FDRThresholds = [1, 1e9, 1e9];
    end
    bc.Map(tobj).NrOfFDRThresholds = size(bc.Map(tobj).FDRThresholds, 1);
    bc.Map(tobj).Name = mapname;
    bc.Map(tobj).BonferroniValue = bfv;
    bc.Map(tobj).UnknownValue = -1;
    bc.Map(tobj).VMPData = aft_SampleBox(mobj{oc}, bb, 1, opts.interp);
end

% put content in new object
bc.NrOfMaps = numel(bc.Map);
bvqxfile_setcont(hfile.L, bc);
